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Abstract. The time evolution of thermal and thermo- 
reactive instabilities of primordial gas in the presence of 
ionizing UV radiation is studied. We obtain conditions 
(depending on density, temperature, and UV radiation in- 
tensity) favorable for the formation of a multi-component 
medium. Nonlinear effects, especially those attributable 
to opacity, can play an important role. In certain pa- 
rameter regimes, dramatic, rapid evolution of ionization 
states away from ionization equilibrium may occur even if 
system control parameters evolve slowly and the system 
starts in or near ionization equilibrium. Mathematically, 
these rapid changes may best be understood as manifesta- 
tions of bifurcations in the solution surface corresponding 
to ionization equilibrium. The astrophysical significance 
of these instabilities is that they constitute a nongravita- 
tional mechanism (thermo-reactive instability) for decou- 
pling isolated spatial regions with moderate initial fluctu- 
ations from their surroundings, even under the influence of 
heating processes due to ionizing UV radiation. In particu- 
lar, isolated, decoupled regions with relatively high neutral 
hydrogen density may result, and these could be associ- 
ated with Lyman alpha clouds of high column density or 
with Lyman limit systems. 
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1. Introduction 

An understanding of processes occurring in diffuse gas in- 
teracting with a photoionizing flux is of direct importance 
in modeling a number of astrophysical phenomena, such 
as quasar absorption clouds (Doroshkevich, Miicket, and 
Miiller, 1990; Ferrara & Giallongo, 1996) and the struc- 
ture of neutral hydrogen clouds in galactic halos (Ferrara 
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& Field, 1994). The effects of heating due to photoion- 
ization in galaxy formation have been considered by Efs- 
tathiou (1992) in the context of mechanisms for suppres- 
sion of dwarf galaxy formation and by Navarro and Stein- 
metz (1996) in their study of the "overcooling" problem. 
These authors found that photoionization alone could not 
provide the heating mechanism required, i.e., supernova 
feedback should also be included. Cosmological simula- 
tions of galaxy formation including star formation and su- 
pernova feedback have been carried out by Yepes, Kates, 
Klypin, & Khokhlov (1996; hereafter YKKK). A reason- 
able hypothesis is that photoheating indirectly influences 
star formation by regulating the conditions for thermal in- 
stability and thus the formation of a multiphase (cloudy) 
medium. 

Since many of the observable phenomena involving 
photoionization require a prediction of statistical prop- 
erties of objects, it is important to integrate the local 
dynamics of diffuse gas (i.e., heating, cooling, ionization, 
formation of neutral hydrogen) into an overall approach 
including hydrodynamics and the evolution of large-scale 
structure. As these processes involve a large dynamical 
range of scales, some of them will inevitably occur below 
the limits of resolution of a numerical simulation. Due to 
the presence of nonlinearities, the effects of fluctuations 
below a typical cell size will require special attention, es- 
pecially when instabilities are involved. 

It has long been appreciated (Field 1965; Defouw 1970; 
Balbus 1986; Fall & Rees; 1985; Ibafiez & Parravano 1982) 
that gas subject to cooling (and perhaps heating) pro- 
cesses in a cosmological setting may become thermally un- 
stable, leading to enhanced cooling and the formation of 
a multiphase medium including cool clouds (Bcgclman & 
McKee, 1990). Cool clouds play an important role in the- 
ories of the interstellar medium (McKee & Ostriker, 1977) 
and galaxy dynamics, due to their effects on the energy 
budget (enhanced cooling), star formation and "supernova 
feedback" (YKKK). 
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Corbelli & Ferrara (1995) have demonstrated the exis- 
tence of so-called "thermo- reactive" instability modes for 
gas containing metals in the presence of ionizing radiation. 
Here we will see that instabilities of thermo-reactive type 
are also possible in gas of primordial composition. In par- 
ticular, our results imply that thcrmo-reactivc instabilities 
could play an important role in the formation of at least 
some population of observed Lyman limit systems. 

The principal goal of this paper is to study the dy- 
namics of fluctuations in a gas of primordial composition 
which may already have been compressed to high ambient 
density (compared to the background density of the uni- 
verse) and heated to high temperature as a consequence of 
large-scale structure formation. The intention is to char- 
acterize instability regimes approximately as a function of 
the ambient temperature and density of the gas, with local 
effects of gravitation excluded. Such a characterization is 
expected to be useful for an understanding of conditions 
for star formation in the context of hydrodynamical nu- 
merical simulations. For this application, linear instability 
analysis about ionization equilibrium solutions does not 
tell the whole story: For one thing, even if present, some 
linear instability modes are too slow to be of importance 
on a dynamical timescale. Moreover, some of the instabil- 
ities which occur involve nonlinearity in an essential way 
and could not be detected by linear analysis. (Because 
opacity is involved, the strength of the nonlinear effects 
increases with the size of the perturbed region; see also 
Ferrara and Field (1994).) 

Apart from its importance for the intended applica- 
tions, the nonlinear system studied in this paper is quite 
interesting from a purely mathematical point of view: As 
seen here [and as previously pointed out for example by 
Petitjean, Bergeron, and Puget (1992)], the equations of 
ionization equilibrium exhibit multivalued solutions in cer- 
tain regimes of parameter space. We propose that these 
solutions are best understood in terms of the theory of 
bifurcations (or catastrophes): Bifurcations are typical in 
the equilibria of nonlinear systems, and the mathemati- 
cal theory (see for example Arnol'd, 1979; Chow & Hale, 
1985) gives strong hints and indications of the general be- 
havior to be expected. Near bifurcations, a system can de- 
part rapidly from equilibrium even if the control parame- 
ters vary slowly. (In the present case, the "slowly-varying" 
control parameters are any two thermodynamic variables, 
say temperature and - on a still longer time scale - pres- 
sure.) Hence, even a qualitative understanding of system 
evolution near bifurcations requires a time - dependent 
treatment of the coupled system of ionization, heating, 
and cooling equations. The importance of simulating the 
dynamical equations for evolution of ionization states was 
previously discussed for example by Cen (1992), and a 
dynamical treatment has been incorporated into the hy- 
drodynamic code used in a number of papers (Cen et al. 
1990,1992,1993,1994). However, for numerical efficiency it 



would be useful to know under what circumstances ion- 
ization equilibrium is in fact a good approximation. 

The organization of this paper is as follows: Section 
H gives the model system of equations describing the dy- 
namics of the three species H/, He/, and He// for a gas 
of primordial composition. Section shows that the so- 
lution manifold of the ionization equilibrium equations 
can exhibit the properties of a mathematical bifurcation 
surface (Golubitsky & Schaeffer, 1985). Of special inter- 
est are the solution trajectories of constant pressure, for 
which one sees that evolution through equilibrium solu- 
tions is not always possible. Section Q considers selected 
solutions of the full time-dependent equations exhibiting 
some typical behavior. Section ^characterizes the regimes 
of temperature and density in which instabilities relevant 
to numerical simulations (i.e., those with sufficiently short 
timescales) are most likely to occur. Section ^ summarizes 
and describes in particular some potentially interesting 
implications for Lyman limit systems. 

2. Equations describing thermal and ionization 
evolution 

Consider first the time-evolution of a sufficiently small re- 
gion such that the density p{t) and gas temperature T{t) 
may be treated as approximately uniform. (This approx- 
imation will be discussed below in connection with the 
optical depth t^.) Self-gravity is neglected. The dynamics 
in this region are assumed to depend on Hubble expansion, 
cooling processes, and heating due to a UV background. 
The UV background is assumed to be described by an 
externally generated, explicitly time-dependent, spatially 
uniform flux J{t). 

We restrict our consideration to a region containing 
hot gas with primordial abundances, i.e. the total he- 
lium number density uhc is 10 % of the total hydrogen 
number density ?i// . The medium is characterized in addi- 
tion by the number densities n.///,n////,n//e/,n//e// and 
nHeiii of its ionization states HI, HII, Hel, Hell, and 
Helll, respectively. One also has n// = n/// -t- nnii and 
nHe = nHei + nHeii + nueiii and the equations of charge 
conservation. 

The evolution of these number densities is determined 
by the equations 

■ = -S,HinHi - ^cHinnirie + nHnnea{HII) (1) 



duH 
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dt 



duHc 
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— (.HelinHell ~ ^e,HeIinHeIine 
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In the above equations, ^i and S,e.i [with i ranging over 
(HI,HeI, Hell)] denote the rates of photo- and coUisional 
ionization, respectively, Ue and n^ are the electronic and 
ionic number densities, and a{i) denotes the recombina- 
tion rate from species i to the next lower ionization state. 
The photo-ionization rate ^i for the j-th species is written 
in the form 



and to construct the total heating rate 



Ux) 



hv 



(4) 



where the Vi are the threshold frequencies of each species, 
J„ is the intensity of the diffuse background UV flux in 
units of ergs s~^ cm~^ Hz""'^ sr~^, and (Ti{v) is the pho- 
toionization cross section, which is approximated using the 
expressions given by Kramers (1923). 

In addition to the primary ionizations described by Eq. 
(0), secondary ionizations could contribute to the num- 
ber density evolution (ShuU, 1979; ShuU & van Steenberg, 
1985; Ferrara & Field 1994). However, there is observa- 
tional evidence for a steepening of the UV spectrum at 
the Hell edge (Jakobsen et al., 1994; Ferrara & Gialongo, 
1996). Hence, the energy of most secondaries will be low 
enough that a large fraction, though not all, of the energy 
will be converted to heat. 

The diffuse flux is diminished by absorption due to 
surrounding gas layers (self-shielding). The optical depth 
at a distance x from the surface is given by 



dx 2^ ni{x')(7i/{i) 



(5) 



where in the summation i ranges over all three species for 
^"11011 < V, over Hel and HI for i^Hd < ^^ < i^HoII and over 
HI for z^Hi <v< t'Hci- 

We model the UV flux in Eq. (Q) as a power law spec- 
trum with spectral index a; it may thus be expressed in 
the form 



Ji,{x) = Jo{v/vll) "exp(-T^(a;)) 



(6) 



where Jo is the intensity at the Lyman limit; we take a ~ \ 
in the calculations that follow. 

For later use, it will be convenient to introduce a 
shielding depth for HI deflned by 



Tff, = -log(eff//(-7rJoGf,)), 



(7) 



where G^j is the parameter calculated by Black (1981) 
for the optically thin case. (Shielding depths for Hel and 
Hell can be deflned similarly.) 

A knowledge of the above quantities allows us to obtain 
the heating rates by photoionization 



r,;(x) 



°° J^{x)h{v -Vi) 

47r ; aJ^jdv 

hu 



(8) 



r = Er. 



(9) 



The cooling rate A is calculated according to Black (1981), 
and Compton cooling is included. 

Now, the flux Ji, appearing in the photoionization rates 
(Q) and the heating rates (||) is a monotonically decreasing 
function of x, which differs from an exponential because 
of the x-dependence of t^^ which itself depends on the 
ni{x). As we will see below, even in a region with mod- 
erate gradients and slowly evolving thermodynamic state 
variables (pressure and temperature) the number densi- 
ties n,j of the ionization states can have large gradients 
and can also change rapidly. Nonetheless, because the x 
- dependence of (Q) and (g) is the same, one may always 
define weighted means rii such that the integral in Eq. M) 
is replaced by 



:^^n^a^{i) 



(10) 



where the summation is as above and Eqs. (^ and (||) re- 
main accurate as written. Of course, these weighted means 
Hi still depend on x, and they are not strictly equal to the 
average number densities. 

Here we identify the true average number densities 
with these weighted values and use them in the evolution 
equations (|^), (2), (3). For the optically thin case, this is 
an excellent approximation, whereas in the optically thick 
case there are some quantitative (but not qualitative) in- 
accuracies as discussed in the Conclusions. 

3. Bifurcation structure in the manifold of solu- 
tions for ionization equilibrium 

One can gain some insight into the time-dependent prob- 
lem by studying the system representing ionization equi- 
librium, which is obtained by setting all time derivatives 
to zero in Eqs. (||), (2) and (3). (For the purposes of 
this discussion, the environment is considered static, i.e., 
constant flux intensity Jq, no cosmological expansion.) 
For a range of number densities uh and temperatures 
T, the complete set of possible solutions for nHi,nHei 
and nueii was identifled and found numerically. Each so- 
lution corresponds to some net heating (or cooling) rate 
r — A oc T = dT/dt, and the locus of equilibrium so- 
lutions deflnes a two-dimensional surface embedded in a 
three-dimensional space with coordinates {tih, T, F — A). 
Fig. ^ shows a plot of this surface in the range 
10^^cm~^ < Uh <'2 X lO^^cm"^ and 1 < r4 < 5, where 
r4 is the temperature in units of lO'^K. The flux inten- 
sity was Jo = J21 10~^^ ergs s~^ cm~^ Hz~^ sr~^, with 
J21 ~ 0.05. This flux density provides a typical example of 
a surface containing a classical catastrophe or bifurcation 
fold which appears to terminate in one or more cusps of 
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the generic type described by Whitney (1955): For temper- 
atures near 1.3 < T4 < 2.1 (the precise range depending 
on the density within the regime 5 x 10""* < Uh < 10^"^) 
there are three solutions for a given nn, T. 

The existence of a bifurcation has profound conse- 
quences. First, as a general property, any system is essen- 
tially nonlinear near a bifurcation, and therefore there is 
no hope of understanding the resulting instabilities from 
a linearized analysis. Secondly, as control parameters of 
a (time-dependent) system initially in a state of equilib- 
rium are slowly varied, departures from equilibrium are 
path dependent: if the path passes through a region with 
no folds, the system may evolve quasistatically through 
a family of nearly ionization equilibrium states. However, 
near a bifurcation even a slow change in control param- 
eters forces the system to undergo such rapid evolution 
that time-dependent solutions would be expected to de- 
part substantially from (ionization) equilibrium. 




Fig. 1. 3-dimensional surface plot of [F — A] as a function of 
uh and T in the range 10~'*cm~'^ < uh < 2 x 10"'^ cm"'^ and 
1 < T4 < 5. The surface is folded and contains at least one 
bifurcation (axes as indicated in the lower panel). 



For a qualitative illustration of what to expect from 
time dependent simulations and for further understand- 



ing of the bifurcations, it is useful to study the behavior 
of T = dT/dt (or equivalently F — A) along curves of con- 
stant pressure P but changing T. The isobaric trajecto- 
ries define a one-parameter family of hyperbolae through 
the two-dimensional (jih, T) parameter space. The time- 
evolution of riH and T for this restricted system corre- 
sponds to moving along one of these hyperbolae either 
in the direction of decreasing T (cooling) or increasing T 
(heating) at the rate given by T. 

Fig. shows a family of plots of F — A vs. T along hy- 
perbolae defined by P=const for varying flux coefRcients 
(from top (1) to bottom (4)) J2i = [0.05, 0.01, 0.005, 0.001]. 
As the flux decreases, the curves have deeper minima. Be- 
ginning on the left, as long as F — A < 0, an isobaric sys- 
tem would move along the curve to the right (decreasing 
temperature). For the low-flux curve (4), the temperature 
drops monotonically until thermal equilibrium is reached 
somewhere below 10^ K. The only other noteworthy fea- 
ture in this temperature range is a higher cooling rate near 
n = 1.7. 



< 




2.5 2.0 1,5 

Temperature T4 (in 10'' K) 
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Fig. 2. Energy loss (F — A) for different UV fluxes at constant 
pressure uhT = 0.95 • 10^ K cm~^. 



The next higher flux (J21 =0.005) curve (3) is triple 
valued between T4 — 1.7 and the cusp near T4 — 1.6. 
The three solutions correspond to three different ioniza- 
tion states solving the equilibrium equations. With respect 
to the full time-dependent equations, the two solutions 
with the highest and lowest values of nni are stable (at- 
tractive) , while the middle solution is unstable (repulsive) . 
Moreover, this regime is likely to be reached by a time- 
dependent system, because dT/dt is negative up to the 
cusp. 

The next higher flux curve (2) also has a triple-valued 
regime, but the curve crosses the T axis (thermal equi- 
librium) near T = 1.4 and thus is not forced to the cusp. 
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2.5 2,0 1.5 

Temperature T"* (in 10^ K) 
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Fig. 3. Energy loss (F — A) for various constant pressures 
and constant UV flux J21 — 0.01. Curve (1): P corresponds to 
T4 = 0.95 at uh = 10"^; Curve (2): P = 5 times higher than 
curve (1) Curve (3): P = 10 times higher than curve (1) Curve 
(4): P = 50 times higher than curve (1) 



Nevertheless, three solutions of the ionization equilibrium 
equations are available in the range where the system is 
still cooling. If some perturbation (such as a small varia- 
tion in the flux or external pressure) occurs, it is possible 
for the system to jump to the lower stable solution, cool, 
and evolve to lower temperatures. 

The upper curve in Fig. includes a triple- valued solu- 
tion within the heating regime (see upper right). A system 
initially at 74 < 1 would be heated until reaching the cusp 
near T4 = 1 and could then jump to the lower branch with 
a lower heating rate. 

The curves in Fig. Q are similar to those of Fig. ||, ex- 
cept that the pressure is varied for constant flux. In the 
lower curve (4) , the solutions are single- valued and do not 
enter the heating regime. The next higher curve (3) has a 
very steep gradient near T4 — 2.0, where it passes rather 
close by the bifurcation. The next curve (2) is analogous 
to curve (3) in Fig. y. Following the upper curve (1) from 
left to right, a system at strictly constant pressure would 
stabilize near T4 = 2 before the multi-valued region is 
reached. However, slight departures from constant pres- 
sure could move the system to the multi-valued region 
and permit a rapid change to the other stable branch. 

The behavior of the curves in Figs. 2 and 3 is quite 
typical of systems with a bifurcation. In particular, they 
illustrate the difficulties likely to arise in a numerical sim- 
ulation. 

For all of the multi- valued curves of Figs. 2 and 3, the 
solution with the largest cooling rate is also associated 
with a drastically higher fraction of neutral hydrogen. In 
particular, at a given pressure it is possible for the system 



to have both an optically thin and an optically thick sta- 
ble state. This possibility gives a hint that there could be 
structures associated with the existence of qualitatively 
different phases in neighboring regions with initially sim- 
ilar (but not identical) control parameters. 

4. Evolution and instability of perturbed regions 
in an ambient hot primordial gas 

In the previous section, we saw that an isobaric system 
may cool in such a way that it is forced out of ionization 
equilibrium. Moreover, we wish to study the evolution of 
a system evolving so quickly that the control parameters 
T and uh need not change simply quasistatically, because 
the timescales associated with cooling and/or expansion 
may be short compared to the relaxation time to ioniza- 
tion equilibrium. Both of these problems require solution 
of the time-dependent equations (0- 3). 




1. 



Fig. 4. Evolution of characteristic parameters as a function 
of the redshift 2 for initial conditions Zi — 4, Th = 20, 
Uh = 0.75 X 10"^ cm"^, J2i{zi) = 0.5. (The flux J2iiz) evolves 
according to Fig. hd.) Starting from the lower left of the figure: 
Two bottom lines: shielding depth thi (definition see Eq. (fm 
for the hot (1; dot-dot-dot-dash) and cold (2; dot-dash) com- 
ponents. Long dashed curve (3): ratio (hot-to-cold) of Jeans 
lengths. Solid curve (4): — logjg(n_ffj/cm~'') for the cold com- 
ponent. Dotted curve (5): Tc in units of T4,. Short dashed line 
(6): Th in units of r4. 



We are interested in the evolution of structures in the 
intergalactic medium and particularly in the formation of 
cold clouds. These effects may occur in a variety of envi- 
ronments, including both gas which is still comoving with 
the Hubble flow and gas in structures that have decou- 
pled from the expansion. In this paper, we derive results 
for the former case (comoving); however, as discussed be- 
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Fig. 5. Evolution as in Fig. |j, except that the density is in- 
creased by a factor of 2. Labeling as in Fig. m. 




3.85 3.80 

Kcdshift z 



Fig. 6. Evolution as in Fig. ft except that the initial flux is 
reduced by a factor of 10 to J21 — 0.05. Labeling as in Fig. 0. 



low, the effect of decoupling the gas from the expansion 
would usually tend to enhance the effects studied here. 

In general, the equations of cooling and ionization are 
coupled to the gas dynamical equations, as well as to grav- 
itation. However, the qualitative character of the insta- 
bilities can be studied without the full generality of gas 
dynamics. The "one-cell" dynamics considered here could 
be generalized to include advection within the context of a 
numerical simulation that does include gas dynamics and 
gravitation (see for example YKKK). 

We consider a gas with primordial composition con- 
sisting of two components: 



— an ambient ("hot") component at temperature Th{t), 
with uniform number density Uh (not to be confused 
with uhi the hydrogen density). The gas is assumed 
as explained above to be comoving with the cosmic 
expansion. 

— an initial perturbation (the "cold" component) con- 
sisting of a slightly compressed region with number 
density ric (assumed spatially constant in the per- 
turbed region) evolving in pressure equilibrium with 
the hot component and therefore at a lower tempera- 
ture Tc = ThTih/nc. For concreteness, we consider the 
cold perturbation region to be a flat slab of thickness 
50 kpc. However, subject to the approximation of Eq. 
([10|), the results should apply to a range of geometries. 

We take a ratio of baryons to dark matter of ils = 0.1 in 
both components and ambient density contrast 5 initially 
with respect to an assumed background number density of 
2x10"'^ cm"3(l-|-Zi)3. This implies n,j is initially (2x10""^ 

In many treatments of thermal instability such as Fall 
& Rees (1985), Field (1965), one is primarily interested in 
whether or not a perturbation can decouple. Here, decou- 
pling is only of interest if it occurs on a time scale short 
compared to tgrav oc iHubbic/\/(l + 5). 

Suppose that the gas is initially in ionization equilib- 
rium. Assuming an initial perturbation of 10% at some 
initial redshift Zi, we study instabilities by following the 
evolution during the time igrav During the evolution we 
monitor parameters characterizing possible decoupling of 
the system including 

— the ratio pd Ph 

— the temperatures Tc and Th 

— the fractions nni of neutral hydrogen 

— the optical depths t^ and r^ 

— the "cooling enhancement" ratio Ec/Eh, where for 
each component E is the energy lost (gained) per unit 
volume 

— the multiplicity of the solutions of ionization equilib- 
rium, if applicable 

— the Jeans lengths for each component 

As disussed by Cen (1992) in a more general context, 
an important technical problem is the presence of sev- 
eral varying time scales, including tgiav , the cooling times 
tcooi, and the various time scales associated with the time- 
dependence of the three species uhi , nnei , and nneii in 
Eqs. (0-3). The complete time dependent equations are 
solved at all times using adaptive timestep control taking 
into account the shortest relevant timescale. However, it 
sometimes occurs that the time-dependent solution is close 
to an (ionization) equilibrium solution. It is then computa- 
tionally efficient to use a different algorithm which finds an 
appropriate equilibrium solution and tests that this solu- 
tion is indeed close to the time-dependent solution sought. 
In this case, the timestep is comparable to the dynamical 
time scale idyn, the minimum of tg^av and icooi- 
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Fig. H shows the evolution with redshift of character- 
istic parameters for initial conditions Zi = 4, T^ = 20, 
Uh = 0.75 X 10~^ cm~^, J21 = 0.5. This density corre- 
sponds to 30 times the background density. For these pa- 
rameters, no thermo-reactive instability occurs during a 



grav 



it 



grav 



10 years), and the cold compo- 



time >> t 

nent does not decouple from the environment: The shield- 
ing depths thi for the hot and cold components coin- 
cide in the figure. Both components remain optically thin 
throughout. The curve (3) giving the ratio of the Jeans 
length of the hot component to that of the cold compo- 
nent attains a value of about 5 during the course of the 
simulation and has flattened out considerably by the end 
of the simulation at z = 1.9. Neither component has be- 
come Jeans unstable within the time scale considered. 

Fig. H is for the same parameters as in Fig. 0, except 
that the density is increased by a factor of 2. This figure 
illustrates how an instability may indeed be present (i.e., 
according to linear perturbation theory) whose timescale 
for decoupling however is longer than the dynamical 
timescale. Clearly, some degree of decoupling has taken 
place: the ratio of Jeans lengths (hot to cold) is 10 and 
still growing, and Th is about 4 times Tc- As before, the 
shielding depths thi for the hot and cold components are 
practically indistinguishable. 

In contrast to the above two figures. Fig. M shows a 
clear case of decoupling. The parameters are as in Fig. S, 
except that the initial flux has been reduced by a factor 
of 10 to J21 = 0.05. Even in this case the optical depths of 
the cold and hot component remain very close until about 
z = 3.77. However, subsequently, the optical depth of the 
hot component barely changes, remaining optically thin as 
in Figs, y and |^, whereas the cold perturbation becomes 
optically thick {uhi increases by more than four orders 
of magnitude) and at the same time rapidly cools and 
thus drastically decreases its Jeans length, until it becomes 
gravitationally unstable. Note that the rapid increase in 
neutral hydrogen density at z = 3.75 is not a numerical 
artifact, but represent rapid evolution as in Figs, g and 
y where the solutions of the static ionization equilibrium 
equations are multi - valued. This evolution provides an 
example of nonlinear thermo-reactive instability and re- 
sults in a strong decoupling of the cold perturbation from 
its environment. 

5. Investigation of instability regimes in the pa- 
rameter space of density and temperature 

As mentioned above, one important goal of this paper is 
to determine what regimes of ambient density and tem- 
perature are favorable for eventual star formation in a 
region containing primordial gas. The arguments given 
here apply to "star formation" as it occurs in the context 
of the McKee - Ostriker (1977) theory of the interstellar 
medium. Now, the formation and evolution of the inter- 
stellar medium are certainly complicated processes involv- 



ing various modes of instability on a range of spatial and 
time scales. However, thermal or thermo-reactive instabil- 
ities are expected to be a necessary prerequisite for such 
a multiphase medium to arise in the first place. Hence, 
a characterization of those regimes of ambient gas tem- 
perature and density likely to produce thermo-reactive or 
thermal instabilities will give information on environments 
favorable for star formation. (Since instability depends on 
UV flux density as well as density and temperature, the 
flux will obviously play an important role in regulating 
star formation.) 

To this end, we have simulated the evolution of "cold" 
and warm components in pressure equilibrium as in Sec- 
tion U under a wide range of initial conditions in number 
density and temperature, and also for a variety of real- 
istic flux and redshift combinations. All simulations be- 
gan with an initial gas density perturbation of 10 % and 
were continued until either the colder component cooled 
to below 5000 K or the local Hubble time igrav defined 
above was exceeded. (Some perturbations which would be 
classified as "unstable" according to linear perturbation 
analysis would nonetheless fail to cool within this time.) 
The dimension of the perturbed region was taken to be 
50/(1 + z) kpc. 

We now present detailed results of several simula- 
tions on a finely-spaced, two-parameter grid (temperature- 
number density) . We first consider the case of initial red- 
shift z=0.1 and flux J21 = 0.01: Fig. |^ is a contour plot of 
the ratio pd Ph evaluated at the end of the simulation. The 
horizontal axis is the initial temperature of the hot com- 
ponent in units of 10* K, and the vertical axis represents 
the initial overdensity of the hot component (with respect 
to the mean gas density). The visible contours represent 
regimes in which the cold component actually did succeed 
in cooling. In the present discussion, we regard a final ra- 
tio Pel Ph > 3 as an indicator of decoupling (instability). 
Note that the density ratio could have been higher than 
its final value during the course of the simulation (see Fig. 
g near z = 3.75). One may observe from the figure that 
the regime of instability as characterized by this indicator 
consists of the region roughly defined by the inequalities 
T4 > 10 (a vertical boundary) and 



(11) 

■^ is the average gas density 



"gasanit > (400 + 25T4)?lgas,0 

where n.gas,o ~ 2 • 10^^ cm^ 
today. 

Fig. ^ is a contour plot of the "enhancement" ratio of 
the cold to the hot component (ratio of total integrated 
heat lost per unit mass during the simulation for the re- 
spective components) for the same simulations as in Fig. 
0. Again, although the final enhancement ratio may be 
lower than the maximum during the course of the simula- 
tion, it still may be regarded as an indicator of instability 
and decoupling. The regime of decoupling as defined by 
this indicator (enhancement greater than 1.5) in Fig. |^ 
coincides quite well with that defined by the final density 
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ratio (Fig. |^). Again, the visible contours are plotted only 
if the perturbed gas actually cooled. 

Since the UV flux changes with redshift, the instabil- 
ity regimes and thus the conditions for development of a 
multicomponent medium will evolve. Hence, one cannot 
necessarily estimate instability conditions in the distant 
past simply by extrapolating the low-redshift results to 
take into account the higher background density. We next 
report the results of higher-redshift simulations performed 
as before on a fine-parameter grid, but for the following 
redshift-flux combinations: 

Zi — 3, J21 = 0.1 (see Fig. ^) and 
Zi = 4, J21 = 0.5 and J21 = 1.0. 

These fluxes are all within the constraints provided by ob- 
servations and large-scale structure simulations (Miicket 
et al., 1996; Bechtold, 1994; Wilhger et al., 1994; Bajt- 
lik, Duncan, & Ostriker, 1988; Lu, Wolfe, & Turnshek, 
1991; Haardt & Madau, 1995). It is quite remarkable that 
for these higher initial redshifts, the left (temperature) 
boundaries of the instability regimes are nearly unchanged 
from the low-redshift case and also hardly differ for the 
two flux values ( J21) considered: the temperature limit is 
still roughly characterized by T4 > 10. 




Temperature T, (in 10' K) 



Fig. 7. Contour plot of the ratio pcj ph evaluated at the end 
of the simulation; initial redshift Zi = 0.1 and initial flux 
J21 = 0.01. Horizontal axis; T^ for the hot component. Ver- 
tical axis: initial overdensity of hot component with respect to 
the background density. Contour labels: pcj ph 



At redshift ^i = 3 for an initial flux J21 = 0.1, one can 
see from Fig. || that the bottom (density) boundary of the 
instability regime can be approximated by 



At the initial redshift Zi = 4 for the lower flux value (0.5), 
the regime of instability is roughly characterized by the 
condition 

ngas,init(^j = 4) > (14+ I.IT4) * rigas.O ' (1 + Zif ■ (13) 

For the higher value (J21 = 1.0), the lower density bound 
is shifted upward by an additional factor of about 1.3. 

In order to investigate the general trend of evolu- 
tion with redshift for the asymptotes described above, 
we have analyzed the instability behavior on a coarser 
grid with respect to {Ti,nH,i ) for a sequence of redshifts 
Zi = 0.1, 1,2,3,4,5 and flux values J2i,i- As a model for 
the redshift dependence of the UV flux, we have used the 
results provided by simulations described in Miicket et 
al. (1996a) which reproduce the observations mentioned 
above quite well (see Fig. O). The results are given in 
Tabic ^. For each initial redshift Zi, the asymptotes (as 
deflned previously) are given approximately by an expres- 
sion of the form 



^gas.init y^i) 



Tq — const, 

(no -I- TO0T4) * rigas.o • (1 + z^f 



(14) 



Zi 


J21 


To 


no 


mo 


0.1 


0.008 


12 


400 


5 


1. 


0.034 


10 


89 


3 


2. 


0.061 


8 


44 


1.75 


3. 


0.074 


10 


25 


1.1 


4. 


0.073 


10 


17 


0.65 


5. 


0.06 


10 


11 


0.5 



t-gas.init 



(z, = 3) > (25 + I.8T4) * ngas,0 • (1 + z^f 



(12) 



Table 1. The parameters To, no, mo for the asymptotes defin- 
ing the regions of instability at given initial redshift Zi and flux 
J2i,i in the n — T plane. 



From these simulations, one may deduce that the pa- 
rameter riQ is most strongly determined by requiring that 
the flux be almost entirely shielded: Since we consider a 
planar slice of flxed comoving size xq = 50 kpc, this re- 
quirement would imply that the critical value no is deter- 
mined by the relation t = XQanHi{l + 2/)^, where Zf is 
the redshift for which nni ~ nongas.o = i^h, i-s., when 
hydrogen becomes neutral. The results of Table || roughly 
flt this relation for a critical optical depth on the order of 
T « 10 — 15 (taking into account typical differences be- 
tween Zi and Zf). The values no{l + Zi)'^ in Table |l| vary 
within 20% around the value 400 (the scatter is probably 
attributable to the coarse computational grid of steps in 
T and tih in the numerical calculations). There is also a 
weak (roughly logarithmic) dependence on the flux, as one 
would expect. 

At higher temperatures (T4 > 10) the terms describ- 
ing coUisional ionization become important. The stage 
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at which the flux is shielded (t « 10) and coohng be- 
comes strong enough is therefore shifted to higher densi- 
ties uh- Although the functional dependence is certainly 
more complicated, the slope mo(l + z)^ can be roughly 
estimated as being proportional to the flux J21 up to red- 
shift z = 4. 

We now return to the simulations discussed above for 
Zi = 0.1 and J21 = 0.01. In contrast to the decoupling 
criteria illustrated in Figs. [^Jl, and ||, in which both cold 
and hot components always lost energy, in the low to mod- 
erate density contrast regime it is possible for decoupling 
to occur in which the hot gas component actually gains 
energy, while the cold perturbation loses energy (cools). 
This effect is seen in the contours of Fig. ^ which form a 
rough band in the figure. Above this band, there is a net 
energy loss in both components; below this band, there 
is a net energy gain in both components. The somewhat 
"speckled" behavior within the band is due to the fact 
that energy may be lost and gained within the course of 
a simulation. The existence of this band could be of great 
interest in a scenario in which the spatial dependence of 
shielding is treated more precisely than here. In this case, 
there could be a persistent core of opaque gas (high neutral 
hydrogen) surrounded by hotter, optically thin gas which 
is heated (See Ferrara & Field, 1994). Above the band, 
one would expect the core to increase gradually in size 
(as in a cooling flow), whereas below the band a pertur- 
bation might be expected to evaporate, especially if heat 
conduction is properly included. 

In order to get an idea of the length scales of the struc- 
tures involved, we apply the results of Cowie & McKee 
(1977) to the parameters considered here and obtain an 
expression for typical values of the critical evaporation 
timescale icvap, which involves the length scale. A critical 
evaporation length of 100 pc may then be estimated from 
the condition iovap > icooi- A larger limit is obtained from 
the condition icvap > ^gravj which leads to a critical size of 
1 kpc. This length scale also follows from a different line 
of argument given by Ostriker & Gnedin (1996). 



Fig. [12|is a contour plot showing the final optical depth 
of the cold perturbation (solid line) and the neutral hy- 
drogen number density ratio of cold to hot components 
Rhi as a function of initial ambient density and temper- 
ature (for initial redshift z=0.1 and flux J21 = 0.01 as be- 
fore). The regions (upper right in Fig. |l2|with both opacity 
T > 1 and Rhi >> 1 coincide with the instability regions 
identified in Fig. [t] on the basis of a moderately elevated 
density contrast. They represent the thermo-reactive in- 
stability as seen earlier for one example of evolution in 
Fig. o. The instability manifests itself as a dramatic rise 
in T and Rhi >> 1. The interpretation will be discussed 
below. 
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Fig. 8. Contour plot of enhancement (see text for definition) 
axes as in Fig. M 




10 20 30 

Temperature T^ (in 10* K) 



Fig. 9. Contour plot of the final value of the ratio pc/ph as in 
Figs. R, but for initial redshift z=3 and initial flux J21 =0.1 



6. Conclusions 

We have obtained regimes of thermal and thermo-reactive 
instability for a range of initial redshifts and fluxes (see 
Sect. 1^ and Table ^). In each case, the boundaries of the 
instability regime as defined by independent criteria (en- 
hancement, density contrast, etc.) were mutually consis- 
tent and only moderately dependent on the instability 
thresholds defined. As seen in Section 0, it is remarkable 
that even in the presence of fiux there is a lower tem- 
perature bound to the instability regime that is almost 
independent of the density as long as the lower bound is 
exceeded. The minimum (initial) number density for in- 
stability to occur corresponds to requiring that the flux 
be almost entirely shielded. 
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Fig. 10. Log-log plot of redshift dependence of the UV flux 

obtained from N-body simulations; here the flux refers to pig_ 12. Simulation parameters z, J21 and axes as in Figs. 
/ — 4:71 J21 H _ hjl. Solid contours: final optical depth of cold perturba- 
tion. Dashed contours: neutral hydrogen number density ratio 
RhI = nHi.c/nHi,h- 



scale structure with galaxy formation (YKKK; Steinmetz, 
1996). In such codes, smoothed estimates of the gas den- 
sity and temperature are produced at every timestep. Con- 
sider for the present discussion a medium with primordial 
abundances; suppose that some estimate of the ultraviolet 
flux is available (see e.g. Miicket et al. (1996)). Now, the 
evolution of thermal and thermo-reactive instabilities in- 
volves the ionization states. These in turn depend in prin- 
ciple not only on the current local values of the temper- 
ature, density, and ultraviolet radiation flux, but also on 
the past history. However, considering the computational 
expense of explicitly simulating the time evolution of all 
relevant ionization species in every cell, it is useful to study 
what can be inferred from a knowledge of the current pres- 
sure, density, and ionizing flux alone. The present paper 
shows that it should generally be possible to estimate the 
instability regimes for any desired redshift/flux combina- 
tions by making use of look-up tables compiled "off-line" 
at each redshift as a function of temperature and pressure. 
The trends found here at least roughly characterize con- 
ditions favorable for formation of a multiphase medium as 
a function of redshift and the ambient temperature and 
density. 

Such a characterization of instability regimes would be 
an extremely useful tool in studying the influence of photo- 
heating on star formation, and hence on galaxy evolution. 
The suppression of thermal instability with increasing UV 
flux also illustrates how an antibiasing mechanism might 
arise if one were to take into account explicitly the emis- 
sion and transport of UV radiation by massive galaxies 
and quasars (Efstathiou, 1992; Haardt & Madau, 1995; 
Ferrara & Giallongo, 1996). 

Referring to Figs, g] and |lj, we identify the thermo- 
reactive instability seen in the evolution of the cold per- 



Fig. 11. Simulations as in Figs. Mandra: regimes in which hot 
gas gained net energy, but cold perturbation lost energy. 



In order to estimate the quality of the approximation 
and averaging procedure in Eq. dlQ), we have performed 
a spatially dependent analysis for selected critical cases. 
This analysis shows that the qualitative picture does not 
change and the error is far below one order of magnitude. 
Moderate quantitative changes would be expected chiefly 
with respect to the characteristic sizes obtained. The ad- 
vantages of the simplifying assumption made at the end of 
Section g are that they allow one to obtain insight into the 
qualitative behavior of optically thick systems and that 
the analysis, which involves an exhaustive search of pa- 
rameter space, can be provided for reasonable computing 
times. 

The above characterization of instability regimes fa- 
vorable for star formation is of particular use in the con- 
text of hydrodynamical numerical simulations of large- 
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turbation as a possible mechanism for at least some Ly- 
man limit systems: During the course of the evolution, 
the optical depth of the perturbation increases dramati- 
cally according to our model. For a region of size 50 kpc, 
column densities ranging from 10^^ to about 10^° cm~^ or 
neutral hydrogen densities of up to 5- 10~^ cm~^ would re- 
sult. (Note that the full hydrogen density Uh also increases 
with respect to its initial value.) 

Realistically, the evolution of the thcrnio-reactive in- 
stability leading to Lyman limit systems will need to be 
modeled more precisely taking properly the spatial varia- 
tion of the optical depth and the onset of Jeans instability 
into account. One would expect that the thermo-reactive 
instability would have a tendency to propagate out from 
an initial core. Once the instability has affected a region 
exceeding the Jeans length, as occurs in Fig. y, gravita- 
tional instability would set in. 

The conditions would thus be favorable for protogalaxy 
formation leading to star formation and hence some degree 
of enrichment of the medium with heavy elements at a 
later stage. Thus, one would expect Lyman limit systems 
that had previously formed from primordial gas according 
to our scenario to contain heavy elements when observed 
now. Indeed, evidence for heavy elements in Lyman limit 
systems has been observed (Petitjean & Bergeron, 1990; 
Petitjean et al., 1994). 

In Fig. |l^ one also finds regimes for which the cold 
component is decoupled {Rhi >> 1) but still optically 
thin. These regions could be of interest in more detailed 
models of the formation of Ly-a clouds. 

Since the main goals of this paper were to characterize 
the energy budget of primordial gas in terms of parameters 
(density, temperature) available from numerical simula- 
tions, the question of bifurcations was not addressed here 
in full mathematical detail. However, in view of the po- 
tential importance of bifurcations for a qualitative under- 
standing of observations (such as Lyman limit systems), 
the general mathematical structure of ionization problems 
in a diffuse gas is worthy of further study in its own right. 
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